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Abstract 


Predetermination  of  background  error  covariance  matrix  B  is  challenging  in  existing  ocean 
data  assimilation  schemes  such  as  the  optimal  interpolation  (01).  An  optimal  spectral 
decomposition  (OSD)  has  been  developed  to  overcome  such  difficulty  without  using  the  B 
matrix.  The  basis  functions  are  eigenvectors  of  the  horizontal  Laplacian  operator,  pre-calculated 
on  the  base  of  ocean  topography,  and  independent  on  any  observational  data  and  background 
fields.  Minimization  of  analysis  error  variance  is  achieved  by  optimal  selection  of  the  spectral 
coefficients.  Optimal  mode  truncation  is  dependent  on  the  observational  data  and  observational 
error  variance  and  determined  using  the  steep-descending  method.  Analytical  2D  fields  of  large 
and  small  mesoscale  eddies  with  white  Gaussian  noises  inside  a  domain  with  4  rigid  and  curved 
boundaries  are  used  to  demonstrate  the  capability  of  the  OSD  method.  The  overall  error 
reduction  using  the  OSD  is  evident  in  comparison  to  the  01  scheme.  Synoptic  monthly  gridded 
world  ocean  temperature,  salinity,  and  absolute  geostrophic  velocity  datasets  produced  with  the 
OSD  method  and  quality  controlled  by  the  NOAA  National  Centers  for  Environmental 
Information  (NCEI)  are  also  presented. 
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1.  Introduction 


In  ocean  data  assimilation  (or  analysis),  the  coordinates  (x,  y,  z )  are  usually  represented 
by  the  position  vector  r  with  grid  points  represented  by  r„,  n  =  1,  2,  N,  and  observational 
locations  represented  by  r("1),  m  =  1,2,  . .  .M.  Here,  N  is  the  total  number  of  the  grid  points,  and  M 
is  the  total  number  of  observational  points.  A  single  or  multiple  variables  c  =  (u,  v,  T,  S,  . . .),  no 
matter  two  or  three  dimensional,  can  be  ordered  by  grid  point  and  by  variable,  forming  a  single 
vector  of  length  NP  with  N  the  total  number  of  grid  points  and  P  the  number  of  variables.  For 
multiple  variables,  non-dimensionalization  is  conducted  before  forming  a  single  vector  c  (Chu  et 
al.  2015)  with  “true”,  analysis,  and  background  fields  (cf,  ca,  c/,)  and  observational  data  (c0)  being 
represented  by  N  and  M  dimensional  vectors, 


(1) 


where  the  superscript  ‘F  means  transpose.  The  innovation  (or  called  the  observational  increment 


d  =  (c„-Hc6), 


(2) 


represents  the  difference  between  the  observational  and  background  data  at  the  observational 
points  r(m).  Here,  H  =[/?,„„]  is  an  M'*N  linear  observation  operator  matrix  converting  the 
background  field  Cb  (at  the  grid  points,  r„)  into  “first  guess  observations”  at  the  observational 
points  r'"'1  (Fig.  1). 

The  analysis  error  (£a)  and  observational  error  (s0)  are  defined  by 


(3a) 


which  are  evaluated  at  the  grid  points.  The  two  errors  are  usually  independent  of  each  other, 


(3b) 


Minimization  of  the  analysis  error  variance 
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E2  =  (*X)  ^min 


(4) 


gives  the  optimal  analysis  field  ca  for  the  “true”  field  ct. 

A  common  practice  in  ocean  data  assimilation  (or  analysis)  is  to  use  a  N*M  weight 
matrix  W=  [wnm]  to  blend  c b  (at  the  grid  points  r„)  with  innovation  d  (at  observational  points 
r°"y)  (Evensen,  2003;  Tang  and  Kleeman  2004;  Chu  et  al.  2004a,  2015;  Galanis  et  al.  2006;  Oke 
et  al.  2008;  Han  et  al.  2013;  Yan  et  al.  2015) 

Ca  =Cb  +  Wd  .  (5) 

Minimization  of  the  analysis  error  variance  with  respect  to  weights, 

SE2  /  dWn,n  =  0  •  (6) 

determines  the  weight  matrix 

W  =  BHr(HBHr  +  R)  1 .  (7) 

Here,  B  is  the  N*N  background  error  covariance  matrix;  R  is  the  M*M  observational  error 
covariance  matrix  and  is  usually  simplified  as  a  product  of  an  observational  error  variance  ( e2o ) 
and  an  identity  matrix  I, 

R  =  e20I.  (8) 

Substitution  of  (7)  into  (5)  leads  to  the  optimal  interpolation  (OI)  equation, 

ca=cb+  BHr(HBHr  +  R)  d  ,  (9) 

which  produces  the  analysis  field  ca  from  the  innovation  d.  The  challenge  for  the  OI  method  is 
the  determination  of  the  background  error  covariance  matrix  B. 

An  alternative  approach  is  to  use  a  spectral  method  with  lateral  boundary  (T)  information 
to  decompose  the  variable  anomaly  at  the  grid  points  [c(r„)  -  cdr,,)]  into  (Chu  et  al.  2015), 
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c 


(10) 


X) -c6(rB)  =  Jjr(rB),  =  A(r„)> 


t=i 


where  {fa}  are  basis  functions;  K  is  the  mode  truncation.  The  eigenvectors  of  the  Laplace 
operator  with  the  same  lateral  boundary  condition  of  (c  -  Cb)  can  be  used  as  the  set  of  the  basis 
functions  {fa}  and  written  in  matrix  (Chu  et  al.  2015) 


®={«U  = 


<4(ri) 

••  4(t) 

Mh) 

<f>2(r2)  .. 

0K  (*2  ) 

i  (Tv) 

02(yN)  .. 

4(*v) 

(11) 


For  a  given  mode  truncation  K,  minimization  of  the  analysis  error  variance  (4)  with  respect  to 
the  spectral  coefficients 


dE2Kldak-  0,  k  =  l,...,K 

gives  the  spectral  ocean  data  assimilation  equation  (Chu  et  al.  2015), 

ca=cb  +  FOr  [OFOr]  1  ®Hrd , 

where  F  is  an  TVx  N (diagonal)  observational  contribution  matrix 


(12) 


(13) 


F  = 


/i  o 

0  f2 

0  0 

0  0 

0  0 

0  0 


0 

0 

0 

0 

0 


0  0 
0  0 
0  0 
f  0 

J  n 

0  \ 
0 


0 
0 
0 
0 

\  0 

0  A. 


M 


(14) 


m= 1 


Here,  the  matrices  ®,  F,  and  H  are  all  given  in  comparison  to  the  OI  equation  (9)  where  the 
background  error  covariance  matrix  B  needs  to  be  determined. 

This  spectral  method  has  been  proven  effective  for  the  ocean  data  analysis.  Chu  el  al. 
(2003a,  b)  named  the  spectral  method  as  the  optimal  spectral  decomposition  (OSD).  With  it, 
several  new  ocean  phenomena  have  been  identified  from  observational  data  such  as  a  bi-modal 
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structure  of  chlorophyll-a  with  winter/spring  (February-March)  and  fall  (September-October) 
blooms  in  the  Black  Sea  (Chu  et  al.  2005a),  fall-winter  recurrence  of  current  reversal  from 
westward  to  eastward  on  the  Texas-Louisiana  continental  shelf  from  the  current-meter,  near¬ 
surface  drifting  buoy  (Chu  et.  al  2005b),  propagation  of  long  Rossby  waves  at  mid-depths 
(around  1000  m)  in  the  tropical  North  Atlantic  from  the  Argo  float  data  (Chu  et  al.  2007),  and 
temporal  and  spatial  variability  of  the  global  upper  ocean  heat  content  (Chu  2011)  from  the  data 
of  the  Global  Temperature  and  Salinity  Profile  Program  (GTSPP,  Sun  et  al.  2009). 

The  spectral  mode  truncation  is  the  key  for  the  success  of  the  OSD  method.  It  acts  as  a 
spatial  low  pass  filter  for  the  fields  to  allow  the  highest  wavenumbers  corresponding  to  the 
highest  spectral  eigenvalues  without  aliasing  due  to  the  information  provided  from  the 
observational  network. 

Questions  arise:  Can  a  simple  and  effective  mode  truncation  method  be  developed  to  take 
into  account  of  model  resolution  (i.e.,  total  number  of  model  grid  points)?  What  are  the  major 
differences  between  OI  and  OSD?  What  is  the  quality  and  uncertainty  of  the  OSD  method?  The 
purpose  of  this  paper  is  to  answer  these  questions.  The  remainder  of  the  paper  is  organized  as 
follows.  Section  2  describes  error  analysis.  Section  3  presents  the  steep-descending  mode 
truncation  method.  Section  4  shows  idealized  “truth”  and  “observational”  fields.  Section  5 
compares  analysis  fields  between  OSD  and  OI.  Section  6  introduces  three  synoptic  monthly 
gridded  world  ocean  temperature,  salinity,  and  absolute  geostrophic  velocity  datasets  produced 
with  the  OSD  method  and  quality  controlled  by  the  NOAA  National  Centers  for  Environmental 
Information  (NCEI).  Conclusions  are  given  in  Section  7.  Appendices  A  and  B  briefly  describe 
several  methods  to  determine  the  H  matrix.  Appendix  C  shows  the  determination  of  basis 
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functions.  Appendix  D  presents  the  Vapnik-Chervonkis  dimension  for  mode  truncation. 
Appendix  E  depicts  a  special  B  matrix  for  this  study. 

2.  Error  Analysis 

Low  mode  truncation  does  not  represent  the  reality  well,  while  high  mode  truncation  may 
contain  too  much  noise.  Let  the  truncated  spectral  representation  s k  in  (10)  at  the  grid  points 
form  an  A-dimensional  vector, 

S£  ~  [•''a  (r,  )■>  A  (L  i’ •  ••’  A  (L )]  •  (15) 

The  M-dimensional  innovation  vector  [see  (2)] 

dr=[rf(rl,V(r(2,),...,<V"1)] 

at  observational  points  can  be  transformed  into  the  grid  points 

M 

’Zkj'"' 

D,=D( r,)  =  “L- - ,  /,=£>„,  (16) 

J  n  m=\ 

where  D( rn)  represents  the  observational  innovation  at  the  grid  points, 

D(rn)  =  Cc(rn)-Cb(rn)-  (17) 

Lrom  Eq(3a),  observations  at  grid  points  are  computed  using  c0(rn)  =  Hrc0(rm)  .  The  original 
background  state,  cb( r«),  keeps  in  the  grid  space.  The  matrix  form  of  (16)  is 

FD  =  Hrd ,  (18) 

where  fn  denotes  contribution  of  all  observational  data  unto  the  grid  point  r„.  The  larger  the  value 
of fn,  the  larger  the  observational  influence  on  that  grid  point  (rn).  D  is  an  A-dimensional  vector 
at  the  grid  points, 

D  T=(Dl,D2,...,DN)  (19) 
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The  analysis  error  (i.e.,  analysis  ca  versus  “truth”  ci)  in  the  spectral  data  assimilation  [see  (10)]  is 
given  by 


£aiTn)  =  Ca^n)~CM 

=  [°a  (r„ )  -  (r„ )]  -  ico  (r„ )  -  (r„ )]  +  ico  (r„ )  - ' c,  (r„ )]  (20) 

=  £K(r„)-D(.rn)  +  £o(rn) 

Here,  (10)  and  (17)  are  used.  The  analysis  error  is  decomposed  into  two  parts 

£a(rn)  =  £K(rn)  +  £o(rn)>  (21) 

with  the  truncation  error  given  by 

^(0  =  ^(0-D<X),  (22a) 

and  the  observational  error  given  by 

£o(rn)  =  CX*„)-CA*n)-  (22b) 


3.  Steep-Descending  Mode  Truncation 


The  Vapnik-Chervonkis  dimension  (Vapnik  1983;  Chu  et  al.  2003a,  2015)  was  used  to 
determine  the  optimal  mode  truncation  Kopt.  As  depicted  in  Appendix  D,  it  depends  only  on  the 
ratio  of  the  total  number  of  observational  points  ( M)  versus  spectral  truncation  (. K)  and  does  not 
depend  on  the  total  number  of  model  grid  points  (TV).  This  method  neglects  observational  error 
and  ignores  the  model  resolution.  In  fact,  the  analysis  error  variance  over  the  whole  domain  is 
given  by 


El  =  ([*>„]  =  [4F£r]  +2 RFC,]  +  RFc,]  ,  Rf£„]  = 


M 

TV  ' ' 


(23) 


where  e]  is  the  observational  error  variance  [see  (8)].  Here,  the  observational  error  is 
assumaed  the  same  at  grid  points  as  at  the  grid  points.  This  is  due  to  the  simplification  of  the 
error  covariance  matrix  R  =  e]  I.  The  Cauchy-Schwarz  inequality  shows  that 
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El  <([£^F£,])  +  2^([£jF£r])1/([£jF£„])  +  ([*>.]) 

=  E2k+  2 Ek  Vm  /  Neo  +  (M  /  N)e2a 


(24) 


The  relative  analysis  error  reduction  at  the  mode-A'  can  be  expressed  by  the  ratio 

Yk 

Both  Ek  and  Ek- i  are  large  for  small  K  (low-mode  truncation),  which  may  lead  to  a  small  value 
of  ]>k.  Both  Ek  and  Ek- i  are  small  for  large  K  (high-mode  truncation),  which  also  leads  to  a  small 
value  of  yK.  An  optimal  truncation  should  be  between  the  low-mode  and  high-mode  truncations 
with  a  larger  value  (over  a  threshold)  of  )>k-  This  procedure  is  illustrated  as  follows.  The  values 
(72,  72,  ...,  yKii)  are  calculated  using  (25)  from  a  large  Kb  (say  250).  The  mean  and  standard 
deviaition  of  7  can  be  computed  as, 


=  In 


E\_x  +  2 Ek_x  /  Ne0  +  Me]  /  N 

El  +  2 Ek  kIM  /  Ne.  +  Me 2  /  N 


K  =  2,3,. 


(25) 


r  = 


1 

Kb~  1 


Kb 

YjYk> 


K= 2 


t(rK-r)2. 


K= 2 


(26) 


Suppose  that  the  relative  error  reductions  (72,  73,  ...,  )>kb)  satisfy  the  Gaussian  distribution.  An 
100(1  -  d)%  upper  one-sided  confidence  bound  on  7  is  given  by 

Yth=r+zas,  (27) 


which  is  used  as  the  threshold  for  the  mode  truncation.  Here,  z  is  the  random  variable  satisfying 
the  Gaussian  distribution  with  zero  mean  and  standard  deviation  of  1 .  If  several  y  values  exceed 
the  threshold,  the  highest  mode 

Kopr  =  max(  K )  (28) 

YK^Yth 

is  selected  for  mode  truncation.  After  the  mode  truncation  Kopt  is  determined,  the  spectral 
coefficnets  (a/c,  k  =  1,  2,  ...,  Kopt)  can  be  calculated,  and  so  as  the  truncation  error  variance 
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3.3.  Multi-Platform  Observations 

Let  observation  be  conducted  by  L  instruments  with  different  e(0l)  deployed  at  r/(m,)  (mi  =  1, 


L 

2,  ..,  Ml ;  /  =  1,2,  ...,  L).  The  total  number  of  observations  is  M  =  '^MI .  The  M-dimensional 

;= i 


observational  vector  is  represented  by 


c!  = 


c0  (ri(1) )» c0  (ri(2) )» •  •  • » cQ  (r/^1 } ),  c0  (r2(1) ),  c0  (r2(2) ),...,  c0  (r^Ml } ), . . . , 

.C0(1,L))»C0(rf))v..,C0( 


The  observational  error  variance  is  given  by 

(<Fe„)  =  M,[<f]2  +  M2[e<2>]2  +...  +  ML[ef'f  . 

The  relative  error  reduction  y k  for  mode  truncation  (25)  is  replaced  by 


Yk= 


El,  +  /Net?  +  /  N 


El  +2EkYjJm,  Wei  +  XM,(e«)2  IN 


1  =  1 


1=1 


,  K  =  2,3,. 


(29) 


(30) 


(31) 


After  the  mode  truncation  is  determined,  the  OSD  equation  (13)  is  used  to  get  the  analysis  field. 

4.  “Truth”,  “Background”,  and  “Observational”  Fields 


Consider  an  artificial  non-dimensional  horizontal  domain  (-19  <  jc  <  19,  - 1 5  <  y  <  15)  with 
the  four  curved  rigid  boundaries  (Fig.  2): 


10 


0.3  cos 


sin 


'x} 

V10y 


=  4  =  - 


—  —  0.2  sin 


(x) 

1-cos 

(y) 

=  rj  =  \ 

v5y 

\  1 

-n  12  (west) 
n  12  (east) 

—n  /  2  (south) 
n  12  (north) 


(32) 


The  domain  is  discretized  with  Ax  =  Ay  =  0.5.  The  total  number  of  the  grid  points  inside  the 
domain  (N)  is  3569.  Fig.  3  shows  the  first  12  basis  functions  {<f>k},  which  are  the  eigenvectors  of 
the  Laplacian  operator  with  the  Dirichlet  boundary  condition,  i.e.,  b\  =  0  in  (C3)  of  Appendix  C. 
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The  first  basis-function  0(xn)  shows  a  one-gyre  structure.  The  second  and  third  basis 


functions  02(xn)  and  03(xn)  show  the  east- west  and  north-south  dual-eddies.  The  fourth  basis- 


function  04(xn)  shows  the  east-west  slanted  dipole-pattern  with  opposite  signs  in  the 
northeastern  region  (positive)  and  the  southwestern  region  (negative).  The  fourth  basis-function 
0,  (xn)  shows  the  tripole-pattem  with  negative  values  in  the  western  and  eastern  regions  and 

positive  values  in  between.  The  higher  order  basis  functions  have  more  complicated  variability 
structures. 

Two  “truth”  fields  for  the  non-dimensional  domain  with  4  rigid  and  curved  boundaries  (Fig. 
2)  contain  multiple  mesoscale  eddies  (treated  as  “truth”)  given  by 


ct  (x,  y)  =  25  -  y 2  /  40  +  3  cos  [LJ(x,  y)]  sin  \_Lyrj(x, y)  +  /?] 

u 


<E  = - 0.3  cos 

10 


y 

v8y 


sin 


vlOy 


,  r)  =  -0.2  sin  ]  1-cosf  — 


(Lx,Ly,j3)  =  (  3, 2,^/2) 
for  the  large-eddy  field  (Fig  4a)  and  given  by 

c,(*,y)  =  25  -y2  / 40  +  3cos[4<^(x,y)]cos[Ly/7(x,y)  +  /?] 

7  =  £-0.2sin(j|l-cos(j 

=  (7,5,0) 

for  the  small-eddy  field  (Fig.  4b).  The  background  field  is  given  by 


e  * 

i  c  =  — 

-0.3  cos 

sin 

10 

UJ 

lioj 

(33) 


(34) 


ci,(x,y)  =  25-y2  /  40 


(35) 


The  “observational”  points  {r(m)}  are  randomly  selected  inside  the  domain  (Fig.  5)  with  the  total 
number  ( M)  of  300.  The  “observational”  points  {r(m)}  are  kept  the  same  for  all  the  sensitivity 
studies.  The  domain  is  discretized  by  Ax  =  Ay  =  0.5  with  total  number  (N)  of  grid  points  of  3,569. 
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Sixteen  sets  of  “observations”  (c0)  are  constructed  from  Figs.4a,  b  using  the  analytical 
values  plus  white  Gaussian  noises  ( e0 )  of  zero  mean  and  various  standard  deviations  (cr)  from  0 
(no  noise)  to  2.0  with  0.1  increment  from  0  to  1.0  and  0.2  increment  from  1.0  to  2.0  (total  16 
sets),  generated  by  the  MATLAB, 

c0(r'"l)  =  c,(rl"l)  +  «<,(r«">).  (36) 

Figs.  6a  and  6b  show  6  out  of  the  16  constructed  sets  with  a  =  (0,  0.2,  0.5,  10.,  1.6,  2.0).  Both 
OSD  and  OI  methods  are  used  to  get  the  analysis  field  c«(rn)  from  these  “observations”.  The 
bilinear  interpolation  (see  Appendix  B)  is  used  for  the  observation  operator  H  in  this  study. 

5.  Comparison  between  OSD  and  OI 
a.  OSD  Analysis  Fields 

The  steep-descending  mode  truncation  Kopt  depends  on  the  user- input  parameter  e0  [see 
(25)]  and  observational  noise  a.  E2a  and  }’k  are  computed  from  the  “observational”  data  in  Figs. 

6a  and  6b.  The  threshold  of  mode  truncation  (27)  varies  with  the  significance  level  a.  In  this 
study,  (eo,  a)  vary  between  0  and  2;  a  has  two  levels  of  (0.05,  0.10)  with  zo.os  =  1.645,  zo.io  = 
1.287  in  (27).  For  given  values  of  eo  (=  0.2)  and  o  (=  0.8),  the  optimal  mode  truncation  depends 
on  the  significance  level  a  with  Kopt=  58  for  a  =  0.05  (Fig.  7a)  and  Kopt=  67  for  a  =  0.10  (Fig. 
7b).  Most  results  shown  in  this  section  is  for  a  =  0.05  since  it  it  a  commonly  used  significance 
level. 

For  the  large-eddy  field,  Kopt  is  not  sensitive  to  the  values  of  a  and  e„.  It  is  7  in  the 
upper-left  portion  and  6  in  the  lower-right  portion  of  Table  1.  For  the  small-eddy  field,  Kopt 
takes  (58,  67)  for  the  most  cases,  178  for  the  high  noise  levels  (<r>1.8)  and  low  e0  values 
( ea  <  1 .0 ),  and  82  for  the  low  noise  levels  ( a  <  0. 1 )  and  low  e0  values  ( eo  <  0.3 )  (Table  2). 
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The  analysis  field  using  the  OSD  data  assimilation  (13)  for  a  particular  user-input  parameter 
e„  and  noise  level  a,  c°aSD ixn,e0,<j) ,  is  represented  in  Fig.  8a  (the  large-eddy  field)  using 

“observations”  in  Fig.  6a  (with  various  a),  and  in  Fig.  8b  (the  small-eddy  field)  using 
“observations”  in  Fig.  6b  (with  various  a).  Comparison  between  Figs.  8a,  b  and  Figs.  4a, b 
demonstrates  the  capability  of  the  OSD  method  with  the  analysis  fields  c°SD(rn,a,e0)  fully 

reconstructed  for  all  occasions. 

b.  OI  Analysis  Fields 

With  the  assumption  that  the  c  field  is  statistically  stationary  and  homogeneous,  the  OI 
equation  (9)  with  the  R  and  B  matrices  represented  by  (8)  and  (El)  [see  Appendix  E]  is  used  to 
analyze  the  “observational”  data  with  three  user-defined  paramters:  (ra,  n,  e0).  Here,  ra  and  n  are 
the  decorrelation  scale  and  zero  crossing  (n  >  ra);  e0  is  the  standard  deviation  of  the 
observational  error.  Let  these  paramters  take  discrete  values  with  total  number  of  Pa  for  ra,  Pb  for 
n,  and  Pe  for  e0.  In  this  study,  we  set  Pa  =  Pb  =  Pe  =  5.  e0  has  5  values  (0.2,  0.5,  1.0,  1.5,  2.0). 
Considering  the  horizontal  domain  from  -15  to  15  in  both  ( x ,  y)  directions,  ra  takes  5  values  (2,  3, 
4,  5,  6);  (n  -  ra)  takes  5  values  (0.5,  1.0,  1.5,  2.0,  2.5).  There  are  125  combinations  of  (ra,  n, 
e0 )  for  the  test. 

The  analysis  field  from  the  OI  data  assimilation  (9),  (rn,<r,ra,rb,e0) ,  with  four  different 

sets  of  user-input  parameters  (ra,  n  -  ra,  e0 ):  (2,  2.5,  1),  (4,  5.5,  1),  (6,  8.5,  1),  and  (6,  8.5,  2),  are 
presented  in  Fig.  9a  (the  large-eddy  field)  using  “observations”  in  Fig.  6a,  and  in  Fig.  9b  (the 
small-eddy  field)  using  “observations”  in  Fig.  6b.  Comparison  between  Figs.  9a,  b  and  Figs. 
4a, b  demonstrates  strong  dependence  of  the  OI  output  on  the  selection  of  the  parameters  (ra,  n, 
e„).  For  the  large-scale  eddies  (Fig.  9a),  the  analysis  fields  ca  are  very  different  from  the  “truth” 
field  ct  for  ra  =  2,  n  =  2.5,  e„  =  1  for  all  “observations”  (Fig.  6a);  the  difference  between  the 
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reconstructed  and  “truth”  fields  decreases  as  ra  and  n  increase;  the  two  fields  are  quite  similar 
when  ra  =  6,  n  =  8.5  for  both  e0  =  1  and  2.  Such  similarity  reduces  with  increasing  e„.  For  the 
small-scale  eddies  (Fig.  9b),  the  analysis  fields  ca  are  totally  different  from  the  “truth”  field  c,  for 
ra  =  6,  n,  =8.5,  e„  =  1  and  2  for  all  “observations”  (Fig.  6b),  less  different  as  ra  and  n,  decrease; 
and  are  quite  similar  to  ct  when  ra  =2,  n  =  2.5,  e„  =  1 . 
c.  Root  Mean  Square  Error 

The  analysis  field  from  OSD,  c°aSD ,  depends  only  on  the  observational  error  variance  e] 
and  its  uncertainty  is  represented  by  the  root  mean  square  error  ROSD, 


D  OSD  /  \ 

R  (o-,e0)  = 


(37a) 


Average  over  all  the  values  of  e0  leads  to  the  overall  uncertainty 


= J^rEthr<vCT.,j-c,<o]2 


n= 1 


(37b) 


The  analysis  field  using  01  ( c°! )  depends  on  three  user-defined  parameters  (ra,  n,  e0).  Its 
uncertainty  due  to  a  particular  parameter  is  represented  by 


ROI(°,ra)  = 

R0\cT,rb)  =  , 


Lpp 

y  iy/rbre  rb  eQ  n= 1 

\wT (r« ><T>ra>rh>e0)-Vt (r  )]' 

Up  ZZZl 

{  ^rare  ra  eQ  n= 1 

Y™  (rn,cj,ra,rb,e0)-vMntf 

(38a) 


(38b) 


ROI(^eo)  =  , 


1 


'  NP  R 


rn  rh  n= 1 


OI 


(r„ ,  cr,  ,  r6,  e0)  -  ^(r„ )]' 


(38c) 


which  are  compared  to  ROSD(cr)  and  ROSD(cr,eo) . 
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Fig.  10  shows  the  comparison  between  R0I(cr,ra)  and  ROSD(cr)  for  5  different  ra  values: 
(2,  3,  4,  5,  6)  and  two  types  (the  large-scale  and  small-scale)  of  the  “observational”  field. 
ROI(cr,ra )  monotonically  increases  with  a  and  is  generally  larger  than  R0SI> (a) .  For  the 

“observations”  representing  the  large-scale  eddy  fields  (Lx  =  2,  Ly  =  3,  see  Fig.  6a),  ROSD(cr ) 
increases  slightly  from  0.32  for  cr  =  0  to  0.34  for  cr  =  2.0.  However,  R0I(cr,ra  =  2)  is  always 
larger  than  7?OSZ)(cr)and  increases  from  0.37  for  a  =  0  to  1.13  for  a  =  2.0;  R01  (cr, r  >  3) is 
smaller  than  Rosn(cr)  for  small  cr,  equals  R0SI> (a)  at  certain  cro,  and  larger  than  ROSD(cr)  for  a  > 
co.  The  value  of  cro  increases  with  ra  from  0.4  for  ra  =  3  to  1.0  for  ra  =  6.  R0/  ( cr,  r  =  6)  increases 
from  0.13  for  a  =  0  to  0.62  for  a  =  2.0.  For  the  “observations”  representing  the  small-scale 
eddy  field  ( Lx  =  5,  Ly  =  7,  see  Fig.  6b),  Rosn(cr)  increases  slightly  from  0.22  for  a  =  0  to  0.27 
for  a  =  0.4;  evidently  from  0.27  for  o  =  0.4  to  0.40  for  a  =  0.5;  and  slowly  from  0.40  for  a  =  0.5 
to  0.71  for  o  =  2.0.  However,  R01  {o,ra )  is  much  larger  than  Rosn(cr)  for  any  ra.  For  example, 

ROI(<j,ra  -  2)  increases  from  0.43  for  cr  =  0  to  1.14  for  a  =  2.0;  ...,  ROI(cr,ra  -6)  increases 
from  0.89  for  cr  =  0  to  1.06  for  cr  =  2.0. 

Fig.  11  shows  the  comparison  between  R°\<j,rb)  and  Rosd((t)  for  5  different  (n,  -  ra) 
values:  (0.5,  1.0,  1.5,  2.0,  2.5)  and  two  types  (large-scale  and  small-scale)  of  the  “observational” 
fields.  R0I(<7,rb)  monotonically  increases  with  cr  and  is  generally  larger  than  ROSD(cr).  For  the 
“observations”  representing  the  large-scale  eddy  fields  (Lx  =  2,  Ly  =  3,  see  Fig.  6a), 
ROI(cr,rb  —  r  )  monotonically  increases  with  c  from  around  0.2  for  cr  =  0  to  around  0.78  for  cr  = 

2.0  for  all  the  values  of  (n  -  ra)  with  cro  from  0.4  for  (n  -  ra)  =  0.5  to  0.6  for  (n-,  -  ra)  =  2.5.  For 
the  “observations”  representing  the  small-scale  eddy  fields  ( Lx  =  5,  Ly  =  7,  see  Fig.  6b), 
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ROI(a,rb  -ra)  is  much  larger  than  ROSD(cr )  for  any  (n,  -ra)  and  cr.  For  example, 
ROI(a,rb  —  ra=  0.5)  increases  from  0.53  for  a  =  0  to  1.00  for  a  =  2.0;  ...,  ROI(<7,rb  —ra=  2.5) 
increases  from  0.58  for  a  =  0  to  1 .00  for  a  =  2.0. 

Fig.  12  shows  the  comparison  between  R01  (cr,eo)  and  R0SD  ( cr,  ea )  for  5  different  e„ 
values:  (0.2,  0.5,  1.0,  1.5,  2.0)  and  two  types  (large-scale  and  small-scale)  of  the  “observational” 
fields.  First,  ROI(cr,eo)  monotonically  increases  with  cr  and  is  evidently  larger  than  ROSD (cr, eo ) 

for  all  a  and  e0.  Second,  dependence  of  ROSD(a,e0)  on  a  is  insensitive  to  the  change  of  e„.  For 
the  “observations”  representing  the  large-scale  eddy  fields  ( Lx  =  2,  Lv  =  3,  see  Fig.  6a), 
ROI(<j,e0)  is  close  to  ROSD(a,e0)  for  a  <  1.2,  and  much  larger  than  ROSD(cr,e0)  for  a  >  1.2  with 
e0  =  0.2  and  0.5;  and  vice  versa  with  e„  =  1.0,  1.5,  and  2.0.  R0' (cr,eo  =  2.0)  increases  slightly 
from  0.98  at  o  =  0  to  1.08  at  a  =  2.0  and  is  almost  twice  of  ROSD(cr,eo )  for  all  a.  For  the 
“observations”  representing  the  small-scale  eddy  fields  ( Lx  =  5 ,  Ly  =  7,  see  Fig.  6b),  Rni  ( cr,  eo ) 
is  also  larger  than  ROSD(<7,eo) .  For  example,  ROI(cr,e0  =  2.0)  increases  slightly  from  1.37  at  cr  = 
0  to  1.42  at  a  =  2.0,  which  is  2-3  times  of  ROSD(cr,e0  =  2.0)  for  a  <  1.0. 

The  overall  performance  between  OI  and  OSD  with  various  noise  levels  (a)  can  be 
estimated  by  the  error  ratio, 


=  ™  *V)- 


Ru1(cj) 


Lppp  X  X  X  Z  [Cf  ^  ^  ^  ’  rb  ’  ^  )  -  C^n  )]2  •  (39) 

^rarbre  ra  r„  e„  n= 1 


Fig.  13  shows  the  dependence  of  k(ct)  (evidently  less  than  1)  on  cr  for  the  two  types  (large-scale 
and  small-scale  eddies)  of  the  “observational”  fields  represented  by  Figs  6a  and  6b  with  two 
different  significance  levels  (a  =  0.05,  0.10)  for  the  threshold  of  mode  truncation  in  the  OSD 
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method  (27).  At  a  =  0.05  (Fig.  13a),  for  the  large-scale  eddy  field,  k(ct)  takes  0.71  at  a  =  0; 
fluctuates  with  a;  and  decreases  to  0.57  at  a  =  2.0.  For  the  small-scale  eddy  field,  k(ct)  increases 
monotonically  with  cr  from  0.43  at  a  =  0  to  0.67  at  a  =  2.0.  At  a  =  0.10  (Fig.  13b),  for  the  large- 
scale  eddy  field,  k(<t)  takes  1.17  at  a  =  0;  decreases  monotonically  with  a  to  0.40  at  cr  =  2.0.  For 
the  small-scale  eddy  field,  k(ct)  increases  monotonically  with  a  from  0.36  at  a  =  0  to  0.70  at  cr  = 
2.0.  It  means  that  the  OSD  performs  better  for  the  test  case.  Integration  of  k(ct)  over  the  whole 
interval  of  the  noise  level  [0,  2.0]  yields 


rx  =  0.05  a-  0.1 


0.76  0.72  large-scale  eddy 

0.51  0.59  small-scale  eddy 


(40) 


which  means  that  the  overall  error  for  the  OSD  is  76%  (51%)  of  the  01  error  for  the  large-scale 
(small-scale)  eddy  field  for  a  =  0.05.  The  overall  performance  of  the  OSD  method  is  relatively 
insensitive  to  the  selection  of  the  significance  level  a. 

The  computational  cost  of  the  OSD  and  OI  methods  is  comparable  in  the  test  cases.  In  the 
OSD  method,  the  steep-descending  method  for  mode  truncation  requires  (a)  the  computation  of 
a  large  number  Kb  in  Eq(26)  of  eigenvectors,  (b)  the  construction  and  solution  of  the  OSD 
equation  (13)  can  be  done  once  for  all.  In  the  OI  method,  however,  the  construction  and  solution 
of  the  OI  equation  (9)  must  be  repeated  each  time  background/observations  changes. 

6.  Synoptic  Monthly  Gridded  Temperature  and  Salinity  Fields 
The  OSD  method  is  used  to  to  produce  the  synoptic  monthly  gridded  (SMG)  temperature 
(: T)  and  salinity  ( S)  datasets  (Chu  and  Fan  2016a,  Chu  et  al.  2016)  from  the  two  world  ocean 
observational  ( T ,  S)  profile  datasets  [the  NOAA  national  Centers  for  Environmental  Information 
(NCEI)  ‘s  World  Ocean  Database  (WOD)  and  the  Global  Temperature  and  Salinity  Profile 
Program  (GTSPP)].  The  synoptic  monthly  gridded  absolute  geostrophic  velocity  dataset  (Chu 
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and  Fan  2016b)  is  also  established  from  the  SMG-WOD  (T,  S)  fields  using  the  P  vector  method 
(Chu  1995,  Chu  and  Wang  2003).  These  datasets  have  been  quality  controlled  by  the  NCEI 
professionals  and  are  openly  downloaded  for  public  use  at 
http://data.nodc.noaa.gov/geoportal/rest/find/document?searchText=svnoptic+monthly+gridded 

&f=searchPage.  The  duration  is  January  1945  to  December  2014  for  the  synoptic  monthly 
gridded  WOD  (T,  S)  and  absolute  geostrophic  velocity  fields  and  January  1990  to  December 
2009  for  the  synoptic  monthly  gridded  GTSPP  (T,  S)  fields. 

7.  Conclusions 

Ocean  spectral  data  assimilation  has  been  developed  on  the  base  of  the  classic  theory  of  the 
generalized  Fourier  series  expansion  such  that  any  ocean  field  can  be  represented  by  a  linear 
combination  of  the  products  of  basis  functions  (or  called  modes)  and  corresponding  spectral 
coefficients.  The  basis  functions  are  the  eigenvectors  of  the  Laplace  operator,  determined  only 
by  the  topography  with  the  same  lateral  boundary  condition  for  the  assimilated  variable  anomaly. 
They  are  pre-calculated  and  independent  on  any  observational  data  and  background  fields.  The 
mode  truncation  K  depends  on  the  observational  data  and  a  user  input  parameter  e2o  (i.e., 

observational  error  variance);  and  is  determined  via  the  steep-descending  method. 

The  OSD  completely  changes  the  common  ocean  data  assimilation  procedures  such  as 
OI,  KF,  and  variational  methods,  where  the  background  error  covariance  matrix  B  needs  to  be 
pre-determined  since  the  weight  matrix  W  is  used.  However,  the  OSD  uses  the  spectral  form  to 
represent  the  observational  innovation  at  the  grid  points  [see  (17)].  Minimization  of  the 
truncation  error  variance  leads  to  the  optimal  selection  of  the  spectral  coefficients.  Thus,  the 
background  error  covariance  matrix  B  vanishes  in  the  OSD  procedure  since  the  weight  matrix 
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W  is  not  used.  It  is  contrast  to  the  existing  01  method,  where  the  B  matrix  is  often  assumed  to 
be  stationary  and  homogeneous  with  user-defined  parameters. 

The  capability  of  the  OSD  method  is  demonstrated  through  its  comparison  to  01  using 
analytical  2D  fields  of  large  and  small  mesoscale  eddies  inside  a  domain  with  4  rigid  and  curved 
boundaries  as  “truth”,  and  addition  to  the  “truth”  of  white  Gaussian  noises  with  zero  mean  and 
standard  deviations  (cr)  varying  from  0  (no  noise)  to  2.0  with  0.1  increment  at  randomly  selected 
locations  used  as  “observations.”  A  simple  covariance  function  (Bretherton,  1976)  was  used  for 
the  01  procedure  with  three  user-defined  parameters  (ra,  n,  e0)  taking  5  possible  values  each. 
The  OSD  uses  the  same  value  of  e0.  The  performance  of  OSD  and  01  is  compared  by  (1) 
patterns  for  each  set  of  125  combinations  of  parameters,  (2)  root  mean  square  errors  for 
varying  parameters,  and  (3)  overall  root  mean  square  erros.  The  results  show  that  the  overall 
error  reduction  using  the  OSD  is  evident,  which  is  76%  (51%)  [72%  (59%)]  for  significance 
level  a  =  0.05  (a  =  0.10)  of  the  OI  error  for  the  large-scale  (small-scale)  eddy  field.  In  context  of 
practical  application,  synoptic  monthly  gridded  world  ocean  temperature,  salinity,  and  absolute 
geostrophic  velocity  datasets  have  been  produced  with  the  OSD  method  and  quality  controlled 
by  the  NO  A  A  National  Centers  for  Environmental  Information  (NCEI). 

Two  issues  need  to  be  addressed  on  the  correlation  matrix.  First,  the  comparison  between 
the  OSD  and  OI  is  at  one  particular  instant  in  time.  The  B  matrix  used  in  the  OI  is  based  only  on 
distance.  Second,  in  the  covariance-matrix  based  methods,  when  the  covariance  matrix  is  fixed 
once  and  for  all,  it  is  well-known  that  the  very  first  data  assimilation  cycle  is  doing  well,  but 
subsequent  cycles  are  less  effective  because  the  remaining  error  has  a  tendency  to  be  orthogonal 
to  the  directions  of  the  covariance  matrix.  In  the  OSD  method,  the  correction  is  based  on 
spectral  functions  (i.e.  basis  functions)  chosen  once-and-for  all.  More  sophisticated,  flow-based 
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covariance  matrix  will  allow  01  to  perform  much  better.  Further  verification  and  validation 
under  real-time  ocean  conditions  are  needed  to  verify  the  quality  of  OSD  in  time  cycles  and  to 
compare  between  OSD  and  OI  methods. 


In  the  two  test  cases  (large  and  small  eddy  fields),  it  is  clear  that  the  optimal  mode 
truncation  Kopt  (around  6  for  the  large  eddy  field  and  around  60  for  the  small  eddy  field)  are 
very  closed  to  the  number  of  eigenvectors  required  to  represent  the  truth  field  (Fig.  4).  This 
shows  the  capability  of  the  steep-descending  mode  truncation.  However,  the  performance  of  the 
method  for  the  truth  field  is  a  mixture  of  large  and  small  scales  in  different  parts  of  the  domains 
needs  to  be  further  investigated. 
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Appendix  A.  Determination  of  H-Matrix  Using  All  Grid  Points 

IDW  interpolation,  using  all  grid  points,  is  one  of  the  most  commonly  used  techniques 
for  interpolation  based  on  the  assumption  that  the  value  of  hmn  in  H-matrix  are  influenced  more 
by  the  nearby  points  and  less  by  the  more  distant  points.  Let 

<C  =  ,/(*<"> -x,)2+(/">-^)2  (Al) 

be  the  distance  between  the  grid  point  {pa,  yj)  and  observational  point  (xim\ y(m)).  The  influence  of 
the  grid  point  x„  on  the  observational  point  \(m)  is  given  by  (Spepard  1968) 

*„,=(<cr/£(<cr  <A2> 

n= 1 

where  q  is  an  arbitrary  positive  real  number  called  the  power  parameter  (typically,  q  =  2). 
Another  form  of  hmn  is  given  by  (Franke  and  Nielson  1991) 
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(A3) 


ft mn  ~  N  9 

n= 1 

where  Dim)  is  the  distance  from  the  observational  point  \{m)  to  the  most  distant  grid  point. 
Eq.(A3)  has  been  found  to  give  better  results  than  (A2)  .  As  a  result,  Cb(x(m),  t),  is  somewhat 
symmetric  about  each  grid  point. 

Appendix  B.  Determination  of  H-Matrix  Using  Neighboring  Grid  Points 

Consider  the  position  vector  x  =  (x,  y )  located  inside  the  grid  cell  (Fig.  B-l), 

x,<x<  x,.„  y;.  <  y  <  yj+1 . 

Mathematically,  the  variable  cb  at  r  (inside  the  grid  cell)  can  be  represented  approximately  by  a 
polynomial, 

ci(r)  =  ZZ^(x-x!)“(j;-3;,)/?  (Bl) 

a= 0  J3= 0 

where  L  =  1  refers  to  the  bilinear  interpolation,  and  L  =  3  leads  to  the  bicubic  interpolation.  For 
the  bilinear  interpolation,  Eq.(Bl)  becomes 

Cb  (r)  =  4)0  +  Ao(x  ~xi)  +  A)i(y  -  Jy)  +  A„(x  -  x,)(y  -  y} )  (B2) 


or  in  matrix  notation, 


cb{r)  =  [\ 


(x— X;)] 


A)o  A\ 

i 

_Ao  Ai_ 

_( y-yj)_ 

(B3) 


Since  a,  at  four  neighboring  grid  points:  cb{xi,  yj),  cb(xi+ 1,  yj),  cb(xm,  yj),  cb(xm,  y/+i)  are  given, 
substitution  of  the  four  values  into  (B2)  leads  to  the  determination  of  the  four  coefficients  Aoo, 
A io,  Aoi,  An.  Using  these  coefficients,  the  bilinear  interpolation  (B2)  becomes 
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C,M)  = 


cb{xM,yj+l)  .  v  ^  cb(xi+l,yj) 


(*m  -xi)(yj+i-yj) 


(x-xj)(y-yi)  + 


(xm  ~xi)(yj+ 1  -y,) 


(x-x/)(>’.+l  -y) 


+  - 


i  (xi  ’  y j+i ) 


(xi+i-x)(j-v.)  +  - 


cb(xny,) 


(x,+1  -x)(yJ+l  -  y ) 


(*/+l  -  4  X^+l  -  Ly  )  (*/+l  -  4  )(^+l  -  y} ) 

Let  the  observational  point  r'"'1  be  located  in  the  grid  cell, 

xi  —  x<m)  <  xi+v  yj  —  y(m)  <  L’y-n  • 

Evaluation  of  Cb  at  the  observational  point  r'"'1  using  (B3)  leads  to 

cb{r(m))  =  Pi”')cb(xi,yj)  +  p^cb(xM,y . p^cjx^y /+1)  +  . p(^  cb(xM,yj+l) 


(B4) 


(B5) 


where  the  proportional  coefficients  {  p(m\p('^),p{'H)^,p{")  ^ }  are  defined  by 


(xM-Xi)(yj+1-yj)  ’  MJ  (xM-xtXyJ+l-yj) 


(m) 

T 

i,j+ 1 


(xM-x™x/a)-yj) 

ixM-xi\yj+\-yj) 


(m) 

r 

i +1,7+1 


C XM  -Xi)(yj+ 1  -^y) 


(B6) 


It  is  noted  that  the  proportionality  coefficients  { pm)  }  depend  solely  on  the 

location  of  the  observational  points  (r(w)),  and 

p{m)  +  p(m)  +  p(m)  +  p(m)  =  1 .  (B7) 

r  U  MJ  ij+ 1  MJ+ 1  V  7 

Setting  L  =3  in  (Bl)  leads  to  the  bicubic  spline  interpolation, 


ciXr)  =  4o  +  4o  (*  “  *<0  +  4)i  O’  “  ) +  4 1  (*  “  Xt  )(y  -  yj ) 

+4>o (x  -  xt  f  +  A02  (y  -  yj  f  +  A30 (x  -  x,  )3  (B8) 

+A2l(x-  x(.  )2  (y  -yj)  +  An  (x  -  xf  )(y  -  )2  +  AQ3  (y  -  y.  f 


or  in  matrix  notation, 
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ci(r)  =  [!  (*-*,-)  ( x-xt f 


(B9) 


which  is  rewritten  by 

c6(r)  =  [l  x 

Ao  4i  y 
Ao  An  j 
A20  a21 

-Ao  ® 

Determination  of  the  ten  coefficients  (Aoo,  Aoi,  Am,  A03,  A10,  An,  An,  A20,  A21,  A30)  requires  not 
only  the  values, 

Ao  =  cb(xi>yj)> 

Ao +AoAx+ Ao  ( Ax)2 + Ao  ( Ax)3  =  cb  (*,-+i > yj )  > 

Ao  +  A  A y  +  4)2  (Ay)2  +  AsAyf  =  cb{xt,yj+l) , 

2  2 

A00  +  Al0  Ax  +  y  +  v4uAxAy  +  v420(Ax)  +  A02(Ay) 

+A30(  Ax)3  +  ^421(Ax)2  Ay  +  ^412Ax(Ay)2  +  ^403(Aj;)3  (Bll) 

=  cb(xM,yj+ j), 

but  also  the  derivatives  at  the  neighboring  grid  points 

Ao  =  dcb (■ xi ,yj)!dx  =  [cb (xi+1  ,yj)~  Cb (xM ,  yj )]  /  2 Ax, 

Al=dcb{xi,y])ldy  =  [cb  (x;. ,  yy+1 )  -  cb  (x,. ,  )]  /  2  Ay , 

4o  +  2(Ax)40  +  3(Ax)240  =  3c6(xi+1,y7.)  /  5x 

=  [c*  (*;+2  >  *; )  -  cb  (*,• ,  y2 )]  /  2  Ax, 


(x-x,.)3] 


Ao  4)1  4)2  ^03 


4)2  43 

42  0 

0  0 
0  0 


Ao 

4i 

^12 
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Ao 

^21 
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3  ~ 
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(BIO) 
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Ao  +  (AvMi  +  (Ay)2  ^i2  =  5c*(w.,+ 01  dx 

=  \-cb(xi+vyJ+i)-  cb(xi-vyj+ !)]  1 2Ax . 


Ai  +  (Ax)Au  +  (A x)2A2I  =  dcb(xM,yj)/  dy 
=  [cb  (xM ,  yj+l )  -  cb  (xM ,  yM  )]  /  2 Ay, 

A0 1  +  2(Ay)A0l  +  3(Ay)2  A03  =  dcb  (xt ,  yj+l )  /  dy 
=  [cA (x,- ,  yJ+ 2 )~cb(xi,  yj )]  /  2 Ay. 


(B12) 


The  solution  of  the  above  set  of  10  linear  algebraic  equations  (Bll)  and  (B12)  leads  to  the 
determination  of  the  ten  coefficients  ( Aoo ,  Aoi,  A02,  A 03,  A 10,  Aw,  An,  A 20,  A21,  A 30).  It  is  noted  that 
values  of  Cb  at  the  10  neighboring  grid  points  (xh  yj),  (x,-+i,  yj),  (xh  v/+i  ),  (x;+ 1,  V/+i  ),  (xm,  yj), 
(xj,  yj.  1),  (x/+2,  yj),  (xm,  V/+i),  (x/+i ,  >7-1),  (x/,  V/+2)  are  used  to  solve  (Bll)  and  (B12).  Following 
(BIO),  interpolation  of  Cb  at  the  10  neighboring  grid  points  on  the  observational  rim)  [=(x("°,  y{m)] 
using  the  bi-cubic  interpolation  is  given  by 
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(B13) 


Thus,  an  equation  similar  to  (B5)  can  be  written  for  evaluating  ch  at  the  observational  point  r('M) 
with  the  known  10  coefficients  (Aoo,  Ao\,  A02,  Ao3,Aw,  An,  An,  A20,  A21,  A30, 


cb  (r(m) )  =  pl”)cb  (x«  >yj)+ p™  c>  (x«+i  >yj)+ p^I  cb  (x;  ’  yJ+i ) 

+pi™l,  cb  (xm  .  yJ+ 1 ) + p™  cb  (xi-\  >yj)+ p™  (xi .  vj- 1 )  +  p{™1  cb  (x;+2  >  yj ) 

+F^+I  (T-1  •  Fy+1 )  +  p%  cb  (T+1 >  Ty-i )  +  P™  c„  (*,- ,  yj+2 ) 


(B14) 
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where  the  10  corresponding  coefficients  {j o(m) ,p(m)  ,p(m)t,p(m) 


Am) 


Am) 


U+ 1 


Z+1J+1 


?0) 

i,+2j 


?(m) 

Hi+i 


Am) 


p{m)2 }  are  analytically  determined  and  depends  solely  on  the  location  of  the  observational  points 
(rtm}),  and 


p(m)  +  p(m)  +  p(m)  +  p(m)  +  p(m)  +  p(m)  +  p(m)  +  p{m)  +  p(m)  +  p(m)  = 

r  i,j  *  i+hj  r  ij+l  r  z+1,7+1  *  i-l,j  *  z'J-1  *  i,+2j  r  i-lj+l  *  i+lj-l  J'i,j+  2 


(B15) 


Since  only  10  neighboring  grid  points  are  used  to  interpolate  at  the  observational  point  rim)  using 
the  bicubic  interpolation,  the  matrix  H  has  only  10  non-zero  values  in  each  row.  However,  it  is 
too  tedious  to  write  it  out. 


Appendix  C.  Basis  Functions 

As  pointed  by  Chu  et  al.  (2015),  three  necessary  conditions  should  be  satisfied  in 
selection  of  basis-functions  { <pk(r)  }:  (i)  satisfaction  of  the  same  homogeneous  boundary 

condition  of  the  assimilated  variable  anomaly,  (ii)  orthonormality,  and  (iii)  independence  on  the 
assimilated  variables.  The  first  necessary  condition  requires  the  same  boundary  condition  for  (c 
-  Cb)  and  the  basis  functions  { <pk } .  The  second  necessary  condition  is  given  by 


JJ^(r)4(rMr  =  4t 

r 

where  dkk’  is  the  Kronecker  delta, 

[0  if***' 
Skk'  [1  if  k  =  K 


(Cl) 


(C2) 


Due  to  their  independence  on  the  assimilated  variable  (the  third  necessary  condition),  the  basis- 
functions  are  available  prior  to  the  data  assimilation. 

The  basis  functions  are  the  eigenvectors  { (j>k }  of  the  Laplacian  operator  with  the  same 
boundary  condition  as  the  variable  anomaly  (c  -  cb). 
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V2A=-44,  |A (r )e* V (j)k  +  h2 (r )(f)k ]  |r  =  0,  k  =  1, (C3) 

Here,  {Xk}  are  the  eigenvalues,  e  is  the  unit  vector  normal  to  the  boundary;  r  denotes  a  moving 
point  along  the  boundary,  and  [  h, (r),  b2(r)  ]  are  parameters  varying  with  r.  The  boundary 

condition  in  (C3)  becomes  the  Dirichlet  boundary  condition  when  b\  =  0,  and  the  Neumann 
boundary  conditions  when  In  =  0.  As  pointed  by  Chu  et  al.  (2015),  different  variable  anomalies 
have  different  [  h,(r),  h2(r)  ].  For  example,  the  temperature,  salinity,  and  velocity  potential 

anomalies  have  bi  =  0  for  the  rigid  boundary  and  hi  =  0  for  the  open  boundary.  However,  the 
anomaly  has  hi  =  0  for  the  rigid  boundary  and  bi  =  0  for  the  open  boundary.  It  is  obvious  that  the 
eigenvectors  { <f)k }  are  orthonormal  and  independent  of  the  assimilated  variables. 

Appendix  D.  Vapnik-Chervonenkis  Dimension  for  Mode  Truncation 

The  Vapnik-Chervonenkis  dimension  (Vapnik  1983;  Chu  et  al.  2003a,  2015)  is  to  seek 
the  optimal  mode  truncation  on  the  base  of  the  first  term  of  the  analysis  error  (23), 

jyun-DDf 

n ={[&**])="'  jV_1 -  (on 

with  the  cost  function 


JK  =Jtr+/j(K,M,a ) 
/u(K,M,a)  =  J , 


[\n(2M  /  K)  +  l]-\n(a  /  M) 
M  IK 


K  =  1,2,. ,.,oo 


(D2) 


Here,  a  (  «  1 )  is  the  significance  level.  J*  is  the  upper  bound  of  Jtr.  For  a  given  M,  Jtr  decreases 
monotonically  with  K;  ju  increases  with  K  if  a  is  given.  The  optimal  mode  truncation  is  through 
the  minimization  of  the  cost  function, 

mm(JK)  =  JKpt.  (D3) 
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This  method  neglects  observational  error  [only  first  term  of  (23)  considered]  and  ignores  the 
model  resolution  (represented  by  the  total  number  of  grid  points  TV).  The  ratio  of  observational 
points  (M)  and  the  spectral  truncation  ( K)  is  the  key  to  determine  the  optimal  mode  truncation 


K, 


opt • 


Appendix  E.  B-Matrix 

The  B  matrix  is  often  established  based  on  the  assumption  of  statistical  stationarity  and 
homogeneity  of  the  reconstructed  field  with  a  simple  covariance  function,  for  example 
Bretherton  et  al.  (1976)  proposed 


B  =  M  ,  b.  = 

L  lJ  J  NxN  lJ 


f  r2  ^ 

1-^ 

V  ^  ) 


exp 


(  2\ 
rJL 
r2  , 

\  “  J 


rij  =  |r;  —  rj |  >  rb>ra. 


(El) 


depending  on  distances  only.  Here,  r,j  is  the  distance  between  the  two  grid  points  r,  and  raf,  ray 
and  n  are  the  decorrelation  scale  and  zero  crossing.  To  conduct  the  OI  data  assimilation,  the 
three  parameters  (e0,  ra,  n)  need  to  be  defined  by  user.  Chu  et  al.  (1997,  2002)  compute  auto¬ 
correlation  functions  from  historical  observational  data  to  fit  the  Gaussian  function  and  get  de- 
correlation  scales  for  the  B  matrix.  Recent  studies  show  that  some  variables  such  as  upper  ocean 
current  speed  does  not  satisfy  the  normal  distribution,  but  the  Weibull  distribution  (Chu  2008, 
2009). 
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Table  1.  Dependence  of  Kopt  on  (cr,  e0)  for  the  large-eddy  field  shown  in  Fig.  6a  with 
significance  level  a  =  0.05. 


Table  2.  Dependence  of  Kopt  on  (a,  e0)  for  the  small-eddy  field  shown  in  Fig.  6b  with 


significance  level  a  =  0.05. 
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Figure  Captions 


Fig.  1 .  Illustration  of  ocean  data  assimilation  with  cb  located  at  the  grid  points,  and  c„  located  at 
the  points  The  ocean  data  assimilation  is  to  convert  the  innovation,  d  =  c0  -  H Cb, 

from  the  observational  points  to  the  grid  points. 

Fig.  2.  Horizontal  non-dimensional  domain  with  four  curved  rigid  boundaries  with  each 
boundary  given  by  Eq.  (32). 

Fig.  3.  Basis  functions  from  cj)i  to  <j)i2  for  the  domain  depicted  by  Eq.(32). 

Fig.  4.  “Truth”  field  ct  taken  as  (a)  the  analytical  function  (33)  with  large-scale  eddy  field  Lx=3, 
Ly  =  2,p  =  nil,  and  (b)  the  analytical  function  (34)  with  small-scale  eddy  field  Lx  =  7,  Ly  =  5,  /?  = 
0. 

Fig.  5.  Randomly  selected  locations  (total:  300)  inside  the  domain  as  “observational”  points. 

Fig.  6a.  “Observational”  data  (c0)  from  Fig.  4a  with  added  white  Gaussian  noises  of  zero  mean 
and  various  standard  deviations:  (a)  0  (i.e.,  no  noise)  (b)  0.2,  (c)  0.5,  (d)  1.0,  (e)  1.6,  and  (f)  2.0. 

Fig.  6b.  “Observational”  data  ( c0 )  from  Fig.  4b  with  added  white  Gaussian  noises  of  zero  mean 
and  various  standard  deviations:  (a)  0  (i.e.,  no  noise)  (b)  0.2,  (c)  0.5,  (d)  1.0,  (e)  1.6,  and  (f)  2.0. 

Fig.  7.  Dependence  of  E2a  and  )>k  on  K  for  the  “observational”  data  for  the  small-scale  eddy  field 

with  a  =  0.8  and  e0  =0.2  at  two  significant  levels  of  (a)  a  =  0.05  (z  0.05  =  1.645)  and  (b)  a  =  0.10 
(z  o.io  =  1.291)  as  the  threshold  of  mode  truncation  [see  Eq.(27)].  The  optimal  mode  truncation  is 
58  for  a  =  0.05  and  67  for  a  =  0.10. 

Fig.  8a.  The  analysis  field  ca  obtained  by  the  spectral  data  assimilation  [see  Eq.(13)]  using  the 
steep-descending  mode  truncation  with  the  significance  level  of  a  =  0.05  from  the  “observations” 
shown  in  Fig.  6a  with  6  noise  (a)  levels  (0,  0.2,  0.5,  1.0,  1.6,  2.0)  and  4  values  of  e0'.  (a)  0.2, 

(i.e.,  no  noise),  (b)  0.5,  (c)  1.0,  and  (d)  2.0. 

Fig.  8b.  The  analysis  field  ca  obtained  by  the  spectral  data  assimilation  [see  Eq.(13)]  using  the 
steep-descending  mode  truncation  with  the  significance  level  of  a  =  0.05  from  the  “observations” 
shown  in  Fig.  6b  with  6  noise  (a)  levels  (0,  0.2,  0.5,  1.0,  1.6,  2.0)  and  4  values  of  e0:  (a)  0.2, 
(i.e.,  no  noise),  (b)  0.5,  (c)  1.0,  and  (d)  2.0. 

Fig.  9a.  The  analysis  field  ca  obtained  by  the  OI  data  assimilation  [see  Eq.(9)]  for  “observations” 
shown  in  Fig.  6a  various  noise  levels  with  various  combinations  of  user-defined  parameters  ( ra , 
n,  e0,)\  (2,  2.5,  1),  (4,  5.5,  1),  (6,  8.5,  1),  and  (6,  8.5,  2). 
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Fig.  9b.  The  analysis  field  ca  obtained  by  the  01  data  assimilation  [see  Eq.(9)]  for  “observations” 
shown  in  Fig.  6b  various  noise  levels  with  various  combinations  of  user-defined  parameters  (ra, 
n,  eQ ,):  (2,  2.5,  1),  (4,  5.5,  1),  (6,  8.5,  1),  and  (6,  8.5,  2). 

Fig.  10.  Comparison  between  R0I(cr,ra )  and  ROSD  (cr)  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  varying  parameter  ra  =  (2,  3,  4,  5,  6)  from  top  to 
bottom  with  the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right  panels  using 
“observations”  in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance  level  of  a  = 
0.05;  and  the  dotted  curves  refer  to  the  OI. 

Fig.  11.  Comparison  between  ROI(<j,rb)  and  ROSD(cr)  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  different  (n  -  ra)  =  (0.5,  1.0,  1.5,  2.0,  2.5)  with 
the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right  panels  using  “observations” 
in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance  level  of  a  =  0.05;  and  the 
dotted  curves  refer  to  the  OI. 

Fig.  12.  Comparison  between  ROI(cr,eo)  and  Rosn(<T,eo)  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  varying  parameter  e0  =  (0.2,  0.5,  1.0,  1.5,  2.0) 
from  top  to  bottom  with  the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right 
panels  using  “observations”  in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance 
level  of  a  =  0.05;  and  the  dotted  curves  refer  to  the  OI. 

Fig.  13.  Dependence  of  the  error  ratio  k  [see  Eq.(39)]  on  a  using  “observations”  in  Fig.  6a 
(represented  by  dots)  and  in  Fig.  6b  (represented  by  *)  with  two  different  significance  levels:  (a) 
a  =  0.05,  and  (b)  a  =  0.10. 

Fig.  Bl.  Interpolation  at  an  observational  point  r(m)  from  four  neighboring  grid  points. 


33 


1 

2 

3 

4 

5 

6 
7 


9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35  705 

36  706 

37  707 

38 

3  q  708 

4  0  709 

41  710 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 
61 
62 

63 

64 

65 


Fig.  1.  Illustration  of  ocean  data  assimilation  with  Cb  located  at  the  grid  points,  and  c0  located  at 
the  points  The  ocean  data  assimilation  is  to  convert  the  innovation,  d  =  c0  -  H Cb, 

from  the  observational  points  to  the  grid  points. 
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15 


Fig.  2.  Horizontal  non-dimensional  domain  with  four  curved  rigid  boundaries  with  each 
boundary  given  by  Eq.  (32). 
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Fig.  3.  Basis  functions  from  c|)i  to  4>i2  for  the  domain  depicted  by  Eq.(32). 
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(a)  Lx=3,  Ly=2 
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Fig.  4.  “Truth”  field  ct  taken  as  (a)  the  analytical  function  (33)  with  large-scale  eddy  field  Lx=3, 
Ly  =  2,  /?  =  7i/2,  and  (b)  the  analytical  function  (34)  with  small-scale  eddy  field  Lx  =  7,  Ly  =  5,  /?  = 
0. 
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300  Random  Location 


Fig.  5.  Randomly  selected  locations  (total:  300)  inside  the  domain  as  “observational”  points. 
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(a):  Observation  Data  without  Noise 


(b):  Observation  Data  with  Noise  (  CT=0.2) 
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Fig.  6a.  “Observational”  data  (c0)  from  Fig.  4a  with  added  white  Gaussian  noises  of  zero  mean 
and  various  standard  deviations:  (a)  0  (i.e.,  no  noise)  (b)  0.2,  (c)  0.5,  (d)  1.0,  (e)  1.6,  and  (f)  2.0. 
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(a):  Observation  Data  without  Noise 


(b):  Observation  Data  with  Noise  (  CT=0.2) 


Fig.  6b.  “Observational”  data  (c0)  from  Fig.  4b  with  added  white  Gaussian  noises  of  zero  mean 
and  various  standard  deviations:  (a)  0  (i.e.,  no  noise)  (b)  0.2,  (c)  0.5,  (d)  1.0,  (e)  1.6,  and  (f)  2.0. 
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Fig.  7.  Dependence  of  E2a  and  )>k  on  K  for  the  “observational”  data  for  the  small-scale  eddy  field 

with  a  =  0.8  and  e0  =0.2  at  two  significant  levels  of  (a)  a  =  0.05  (z  0.05  =  1.645)  and  (b)  a  =  0.10 
(z  o.io  =  1.291)  as  the  threshold  of  mode  truncation  [see  Eq.(27)].  The  optimal  mode  truncation  is 
58  for  a  =  0.05  and  67  for  a  =  0.10. 
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Fig.  8a.  The  analysis  field  ca  obtained  by  the  spectral  data  assimilation  [see  Eq.(13)]  using  the 
steep-descending  mode  truncation  with  the  significance  level  of  a  =  0.05  from  the  “observations” 
shown  in  Fig.  6a  with  6  noise  (a)  levels  (0,  0.2,  0.5,  1.0,  1.6,  2.0)  and  4  values  of  e0\  (a)  0.2, 
(i.e.,  no  noise),  (b)  0.5,  (c)  1.0,  and  (d)  2.0. 
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Fig.  8b.  The  analysis  field  ca  obtained  by  the  spectral  data  assimilation  [see  Eq.(13)]  using  the 
steep-descending  mode  truncation  with  the  significance  level  of  a  =  0.05  from  the  “observations” 
shown  in  Fig.  6b  with  6  noise  (a)  levels  (0,  0.2,  0.5,  1.0,  1.6,  2.0)  and  4  values  of  e0\  (a)  0.2, 
(i.e.,  no  noise),  (b)  0.5,  (c)  1.0,  and  (d)  2.0. 
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Fig.  9a.  The  analysis  field  ca  obtained  by  the  01  data  assimilation  [see  Eq.(9)]  for  “observations” 
shown  in  Fig.  6a  various  noise  levels  with  various  combinations  of  user-defined  parameters  (ra, 
n,  e09):  (2,  2.5,  1),  (4,  5.5,  1),  (6,  8.5,  1),  and  (6,  8.5,  2). 
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Fig.  9b.  The  analysis  field  ca  obtained  by  the  01  data  assimilation  [see  Eq.(9)]  for 
“observations”  shown  in  Fig.  6b  various  noise  levels  with  various  combinations  of  user- 
defined  parameters  (ra,  n,  eQ ,):  (2,  2.5,  1),  (4,  5.5,  1),  (6,  8.5,  1),  and  (6,  8.5,  2). 
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Fig.  10.  Comparison  between  7?0/(cr,r )  and  ROSD(cr )  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  varying  parameter  ra  =  (2,  3,  4,  5,  6)  from  top  to 
bottom  with  the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right  panels  using 
“observations”  in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance  level  of  a  = 
0.05;  and  the  dotted  curves  refer  to  the  OI. 
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Fig.  11.  Comparison  between  R0I(cr,rb )  and  R0SD(cr)  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  different  ( n  -  ra)  =  (0.5,  1.0,  1.5,  2.0,  2.5)  with 
the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right  panels  using  “observations” 
in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance  level  of  a  =  0.05;  and  the 
dotted  curves  refer  to  the  OI. 
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Fig.  12.  Comparison  between  ROI(<j,e0 )  and  ROSD(a,e0)  of  the  analysis  fields  from  the  same 

“observations”  with  different  noise  levels  with  varying  parameter  e0  =  (0.2,  0.5,  1.0,  1.5,  2.0) 
from  top  to  bottom  with  the  left  panels  using  “observations”  shown  in  Fig.  6a  and  the  right 
panels  using  “observations”  in  Fig.  6b.  The  solid  curves  represent  the  OSD  with  the  significance 
level  of  a  =  0.05;  and  the  dotted  curves  refer  to  the  OI. 
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